#########################
#   John Henderson      #
#                       #        
#    Amplify:  	        #
#   amplifyUnobserved.R #
#                       #
#   This version:       #
#    Apr 14, 2011       #
#########################

# amplifyUnobserved

# The function 'amplifyUnobserved' implements an traditional sensitivity analysis for matched and unmatched groups with binary treatments 
#  and non-binary (continuous or categorical) responses, given imbalance on unobserved, covariates.  This sensitivity analysis 
#  is based on a one parameter model developed in Rosenbaum (2002). In such a sensitivity test, researchers 
#  assess how robust a finding is to the presence of an unobserved factor that is assumed to be associated with both a treatment 
#  and an outcome, by seeing how that finding changes at varying levels of the confounding association. 
#
# In the model, Gamma denotes the maximum odds-ratio of receiving treatment between any treated unit and any control unit across all 
#  strata. To assess confounding due to a hypothetical unobserved factor, increasing hypothetical odds-ratios are tested successively 
#  until a finding is no longer robust. However, in observational studies, both observed and unobserved factors frequently remain 
#  a problem in the analysis. 'amplifyUnobserved' allows research to see how robust their findings are to remaining confounding due 
#  due to unobserved covariates. The function uses Wilcoxons rank sum statistic, and allows the researcher to report either probability 
#  values or Hodges-Lehmann point estimates of an additive treatment effect (See Rosenbaum (2002) for more details).
#
#  The major difference ihere is that that function loops over fractions and stops at the HL value or propability value no longer 
#   indicates a different from zero result. 
#
# 'amplifyUnobserved' takes five arguments: 
#   treated:  a vector of non-binary responses for treated units 
#   control:  a vector of non-binary responses for control units 
#   Xmat_tr:  a vector or matrix of observed covariates for treated units 
#   Xmat_ct:  a vector or matrix of observed covariates for control units 
#   level:  a vector of quantiles in the distribution of the odds-ratios across treated and control units to be used to calculate bounds
#   matched:  a toogle to indicate whether the data is pair-matched or unmatched 
#   pval:  a toggle to indicate whether the function reports probability values based on a Wilcoxon rank sum test, or 
#          Hodges-Lehmann point estimates of an additive effect 
#   alternative:  this indicates the area in which the hypothesis test will be conducted; can be either "greater" or "less" 
#   sgn:  this indicates that the control vector is greater than the treated vector outcomes, and thus are reversed; bad, but quick solution
#
# 'amplifyObserved' reports a table of probability values or Hodges-Lehmann (HL) point estimates for a treatment 
#  effect at each level of Gamma and Delta.  Note that the HL point estimates take quite a bit longer to compute.  


amplifyUnobserved=function(treated=NULL,control=NULL,matched=T,sig=.1,eps=.001,inc=.0001,digits=5,HLstop=F,...){
 	
	alternative="greater"
	
	tr=treated
	ct=control
				
	# to ensure that the proper distribution is used to estimate Gamma
	if(mean(treated) >= mean(control)){
		alternative='greater'
		sgn=1
	} else if(mean(treated) < mean(control)){
		alternative='greater'
		control=tr
		treated=ct
		sgn=-1
	}
	
	Delta=10e100
	

	if(HLstop==T){
		len=0
#Gamma=c(-9,-5,0)#; starting region
		
#while(len==0){
#Gamma=Gamma+10
#if(matched==T){	
#out=amplifySignrank(treated=treated,control=control,Gamma=Gamma,Delta=Delta,amplify=T,inc=inc,digits=digits,alternative=alternative)
#} else if(matched==F){
#out=amplifyRanksum(treated=treated,control=control,Gamma=Gamma,Delta=Delta,digits=digits,alternative=alternative)
#}
#len=length(which(sign(out$HL*sgn)!=sgn))
#}
		
#in0=Gamma[min(which(sign(out$HL*sgn)!=sgn))-1]
#in1=Gamma[min(which(sign(out$HL*sgn)!=sgn))]
		
#if(in1!=1){
#if(length(in0)==0){
#in0=in1-5
#}
#} else if(in1==1){
#in0=in1
#}		
		in0=1
		in1=100
		while(in1-in0>eps){
			Gamma=in0+(in1-in0)/2
			if(matched==T){	
				out=amplifySignrank(treated=treated,control=control,Gamma=Gamma,Delta=Delta,amplify=T,inc=inc,digits=digits,alternative=alternative)
			} else if(matched==F){
				out=amplifyRanksum(treated=treated,control=control,Gamma=Gamma,Delta=Delta,digits=digits,alternative=alternative)
			}
			if(sign(out$HL)==sgn){
				in0=Gamma
				in1=in1
			}else if(sign(out$HL)!=sgn){
				in0=in0
				in1=Gamma
			}
		}
		
		Gamma=mean(in1,in0)
		if(matched==T){	
			out=amplifySignrank(treated=treated,control=control,Gamma=Gamma,Delta=Delta,amplify=T,inc=inc,digits=digits,alternative=alternative)
		} else if(matched==F){
			out=amplifyRanksum(treated=treated,control=control,Gamma=Gamma,Delta=Delta,digits=digits,alternative=alternative)
		}
	} else if(HLstop==F){
#len=0
#Gamma=c(-9,-5,0)#; starting region
		
#while(len==0){
#Gamma=Gamma+10
#if(matched==T){	
#out=amplifySignrank(treated=treated,control=control,Gamma=Gamma,Delta=Delta,amplify=T,inc=inc,digits=digits,alternative=alternative)
#} else if(matched==F){
#out=amplifyRanksum(treated=treated,control=control,Gamma=Gamma,Delta=Delta,digits=digits,alternative=alternative)
#}
#len=length(which(out$pvalues$one_tailed>sig))
#}

#in0=Gamma[min(which(out$pvalues$one_tailed>sig))-1]
#in1=Gamma[min(which(out$pvalues$one_tailed>sig))]
		
#if(in1!=1){
#if(length(in0)==0){
#in0=in1-5
#}
#} else if(in1==1){
#in0=in1
#}
#pval=1
#hl_sign=-sgn
		in0=1
		in1=100
		while(in1-in0>eps){
			Gamma=in0+(in1-in0)/2
			if(matched==T){	
				out=amplifySignrank(treated=treated,control=control,Gamma=Gamma,Delta=Delta,amplify=T,inc=inc,digits=digits,alternative=alternative)
			} else if(matched==F){
				out=amplifyRanksum(treated=treated,control=control,Gamma=Gamma,Delta=Delta,digits=digits,alternative=alternative)
			}
			if(out$pvalues$one_tailed<sig){
				in0=Gamma
				in1=in1
			}else if(out$pvalues$one_tailed>=sig){
				in0=in0
				in1=Gamma
			}
#hl_sign=out$HL
#pval=out$pvalues$two_tailed
		}
		
		Gamma=mean(in1,in0)
		if(matched==T){	
			out=amplifySignrank(treated=treated,control=control,Gamma=Gamma,Delta=Delta,amplify=T,inc=inc,digits=digits,alternative=alternative)
		} else if(matched==F){
			out=amplifyRanksum(treated=treated,control=control,Gamma=Gamma,Delta=Delta,digits=digits,alternative=alternative)
		}
		
	}
		
	out$HL=sgn*out$HL
	return((out))
	
} 

#END amplifyObserved


#######################
#   John Henderson    #
#                     #
#    Amplify:	      #
#   amplifySignrank.R #
#                     #
#   This version:     #
#    Mar 4, 2011      #
#######################

# amplifySignrank

# The function 'amplifySignrank' implements a two paramter sensitivity analysis for matched groups with binary treatments 
#  and non-binary (continuous or categorical) responses.  This sensitivity analysis is based on Gastwirth et al (1998), 
#  which establishes bounds for matched pairs for binary and non-binary outcomes.  (For non-binary outcomes, Gastwirth bounds 
#  are NOT identcal with the bounds identified in an amplification approach developed in Rosenbaum and Silber (2009) 
#  and Rosenbaum (2002).  Thus, both are allowed here.)  
#
# 'amplifySignrank' takes eight arguments: 
#   treated:  a vector of non-binary responses for treated units (length of treated must equal length of control)
#   control:  a vector of non-binary responses for control units (length of treated must equal length of control)
#   mdata:  a matched data object produced by the 'Match' function in the 'Matching' package prouced by Jas Sekhon
#           (see 'Matching')
#   Gamma:  vector of sensitivity values to check for the association between U and treatment 
#   Delta:  vector of sensitivity values to check for the association between U and response 
#   amplify:  indicator that toggles whether or not amplification or Gastwirth bounds are reported 
#   pval:  indicator that toggles whether probability values or Hodges-Lehmann point estimates are reported 
#   alternative:  this indicates the area in which the hypothesis test will be conducted; can be either "greater" or "less" 

# 'amplifySignrank' reports a table of probability values or Hodges-Lehmann point estimates for a treatment effect at 
#  each level of Gamma and Delta


amplifySignrank=function(treated=NULL,control=NULL,Gamma=c(1:5),Delta=c(1:5),
    amplify=TRUE,inc=.0001,digits=5,alternative="greater",...){

		hlSolve=function(R,Z,t=NULL,matched=F,inc=.0001){
			
			m=length(which(Z==1))
			n=length(R)
			invisible(library(stringr))  
			dg=abs(str_length((inc))-3)
			
			#matched == F is rank sum statistic 
			#matched == T is sum sign rank statistic 
			
			if(matched==F & is.null(t)){
				t=m*(n+1)/2
			}
			if(matched==T & is.null(t)){
				t=m*(n/2+1)/4
			}
			
			p_tau=m_tau=0
			np0=np1=nm0=nm1=1
			
			while(sign(np0)==sign(np1) & sign(nm0)==sign(nm1)){
				if(matched==F){
					np0=sum(rank((R-p_tau*Z))[which(Z==1)]) - t
					p_tau=p_tau+inc
					np1=sum(rank((R-p_tau*Z))[which(Z==1)]) - t
					
					nm0=sum(rank((R-m_tau*Z))[which(Z==1)]) - t
					m_tau=m_tau-inc
					nm1=sum(rank((R-m_tau*Z))[which(Z==1)]) - t
					
				} else if(matched==T){
					np0=sum(rank(abs(((R-p_tau*Z)[Z==1])-((R-p_tau*Z)[Z==0])))[which(((R-p_tau*Z)[Z==1])>((R-p_tau*Z)[Z==0]))]) - t
					p_tau=p_tau+inc
					np1=sum(rank(abs(((R-p_tau*Z)[Z==1])-((R-p_tau*Z)[Z==0])))[which(((R-p_tau*Z)[Z==1])>((R-p_tau*Z)[Z==0]))]) - t
					
					nm0=sum(rank(abs(((R-m_tau*Z)[Z==1])-((R-m_tau*Z)[Z==0])))[which(((R-m_tau*Z)[Z==1])>((R-m_tau*Z)[Z==0]))]) - t
					m_tau=m_tau-inc
					nm1=sum(rank(abs(((R-m_tau*Z)[Z==1])-((R-m_tau*Z)[Z==0])))[which(((R-m_tau*Z)[Z==1])>((R-m_tau*Z)[Z==0]))]) - t
				}
			}
			
			if(abs(nm1) + abs(nm0) < abs(np1) + abs(np0)){
				tau=-median(c(-c(m_tau,m_tau+inc),c(p_tau-inc,p_tau)))
			} else if(abs(nm1) + abs(nm0) > abs(np1) + abs(np0)){
				tau=median(c(-c(m_tau,m_tau+inc),c(p_tau-inc,p_tau)))
			} else{
				tau=0
			}
			
			return(tau)	
		}	
		
	if(is.null(treated) | is.null(control)){
		stop("Need input data")
	}
	if(length(treated)!=length(control)){
		stop("Either this is an unmatched object or this is not 1 to 1 matching") 
	}
	if(length(unique(c(treated,control)))==2){ 
		stop("Responses are binary")
	}
		
	alternative="greater"
	
	if(alternative=="greater" & mean(treated)>=mean(control)){
		tr=treated
		ct=control
		sgn=1
	} else if(alternative=="greater" & mean(treated)<mean(control)){
		tr=treated
		ct=control
		sgn=-1
	} else if(alternative=="less" & mean(treated)>=mean(control)){
		tr=control
		ct=treated
		sgn=1
	} else if(alternative=="less" & mean(treated)<mean(control)){
		tr=control
		ct=treated
		sgn=-1
	}
	
	treated=tr
	control=ct
	
	outs=list()
	outs[[1]]=outs[[2]]=outs[[3]]=lambda=matrix(NA,length(Gamma),length(Delta))
	rownames(outs[[1]])=rownames(outs[[2]])=rownames(outs[[3]])=c(as.character(round(Gamma,digits=3)))
	colnames(outs[[1]])=colnames(outs[[2]])=colnames(outs[[3]])=c(as.character(round(Delta,digits=3)))	
		
	if(mean(tr)>=mean(ct)){
		treated=tr
		control=ct
	} else if(mean(tr)<mean(ct)){
		control=tr
		treated=ct
	}
	
	difs=treated-control
	n=length(difs)
	qs=rank((difs))
	q=rank(abs(difs))
	Tstat=sum(q[which(difs>0)])   
	
	et_hl=array(0,length(Gamma)*length(Delta))
	
	for(i in 1:length(Gamma)){
		
		row.vec=array(NA,length(Delta))
		for(j in 1:length(Delta)){
				
			if(amplify==T){
				lambda[i,j]=(Delta[j]*Gamma[i]+1)/(Delta[j]+Gamma[i])
				gam=log(lambda[i,j])
				del=log(1e100)
					
				prob=lambda[i,j]/(1+lambda[i,j])
				et=sum(q*prob)
				vt=sum(q*q*prob*(1-prob))
				et_hl[(i-1)*length(Delta)+j]=sum(q*prob)
				
			} else if(amplify==F){
				gam=log(Gamma[i]) # 1:Gam
				del=log(Delta[j]) # 1:Del
				
				pie=1/(1+exp(-abs(gam)*(1-0)))
				theta=1/(1+exp(-abs(del)*(2*qs/n)))
					
				if(sign(gam*del)>=0){
					prob=pie*theta+(1-pie)*(1-theta)
				} else if(sign(gam*del)<0){
					prob=pie*(1-theta)+(1-pie)*theta
				}
					
				et=sum(qs*prob)
				vt=sum(qs^2*prob*(1-prob))
				et_hl[(i-1)*length(Delta)+j]=sum(qs*prob)
			}
				
			deviate=(Tstat-et)/sqrt(vt)
			row.vec[j]=1-pnorm(deviate)
		}
			
		outs[[2]][i,c(1:ncol(outs[[2]]))]=c(row.vec)
	}
	
	et_u=unique(et_hl)
	tau_u=array(0,length(et_u))
	#tau_u[1]=round(wilcox.test(treated,control,paired=TRUE,conf.int=TRUE,exact=FALSE)$estimate,digits=digits)
	
	r=c(treated,control)
	z=c(rep(1,length(treated)),rep(0,length(control)))
	
	for(l in 1:length(et_u)){		
		tau_u[l]=hlSolve(R=r,Z=z,matched=T,t=et_u[l],inc=inc)
	}
		
	tau=array(0,length(et_hl))
	
	for(l in 1:length(et_u)){
		tau[which(et_hl==et_u[l])]=tau_u[l]
	}
	
	for(i in 1:length(Gamma)){
		for(j in 1:length(Delta)){
			outs[[1]][i,j]=tau[length(Delta)*(i-1)+j]*sgn
		}	
	}
	
	outs[[3]]=outs[[2]]
	for(j in 1:ncol(outs[[3]])){
		inds1=which(outs[[3]][,j]>.5)
		inds2=which(outs[[3]][,j]<=.5)
		outs[[3]][inds1,j]=2*(1-outs[[3]][inds1,j])		
		outs[[3]][inds2,j]=2*(outs[[3]][inds2,j])
	}

	return(list(HL=round(outs[[1]],digits=digits),pvalues=list(one_tailed=round(outs[[2]],digits=digits),two_tailed=round(outs[[3]],digits=digits)),Gamma=Gamma,Delta=Delta))
} 

#END amplifySignrank


#######################
#   John Henderson    #
#                     #        
#    Amplify:	      #
#   amplifyRanksum.R  #
#                     #
#   This version:     #
#    Mar 4, 2011      #
#######################

# amplifyRanksum

# The function 'amplifyRanksum' implements a two paramter sensitivity analysis for unmatched groups with binary treatments 
#  and non-binary (continuous or categorical) responses.  This sensitivity analysis is based on an amplification approach 
#  of the one parameter model developed in Rosenbaum and Silber (2009) and Rosenbaum (2002).  The function uses Wilcoxons 
#  rank sum statistic, and allows the researcher to report either probability values or Hodges-Lehmann point estimates of an 
#  additive treatment effect (See Rosenbaum (2002) for more details).
#
# 'amplifyRanksum' takes six arguments: 
#   treated:  a vector of non-binary responses for treated units 
#   control:  a vector of non-binary responses for control units 
#   Gamma:  vector of sensitivity values to check for the association between U and treatment 
#   Delta:  vector of sensitivity values to check for the association between U and response 
#   pval:  a toggle to indicate whether the function reports probability values based on a Wilcoxon rank sum test, or 
#          Hodges-Lehmann point estimates of an additive effect 
#   alternative:  this indicates the area in which the hypothesis test will be conducted; can be either "greater" or "less" 
#
# 'amplifyRanksum' reports a table of (approximate) probability values or Hodges-Lehmann (HL) point estimates for a treatment 
#  effect at each level of Gamma and Delta.  Note that the HL point estimates take quite a bit longer to compute.  


amplifyRanksum=function(treated=NULL,control=NULL,Gamma=c(1:5),Delta=c(1:5),
	inc=.0001,digits=5,alternative="greater",...){
	
	if(is.null(treated) | is.null(control)){
		stop("Need input data")
	}
	if(length(unique(c(treated,control)))==2){ 
		stop("Responses are binary")
	}

	# required packages: plyr
	if(length(which(.packages(all.available=T)=='plyr'))==0){
		install.packages("plyr")
	}
	library(plyr)

	alternative="greater"

	if(alternative=="greater" & mean(treated)>=mean(control)){
		tr=treated
		ct=control
		sgn=1
	} else if(alternative=="greater" & mean(treated)<mean(control)){
		tr=treated
		ct=control
		sgn=-1
	} else if(alternative=="less" & mean(treated)>=mean(control)){
		tr=control
		ct=treated
		sgn=-1
	} else if(alternative=="less" & mean(treated)<mean(control)){
		tr=control
		ct=treated
		sgn=1
	}
	
	treated=tr
	control=ct	
	
	q=rank(c(treated,control))
	t=sum(q[1:length(treated)])
	m=length(treated)
	n=length(q)
	qsort=sort(q,decreasing=TRUE)
	
	R=c(treated,control)
	Z=c(rep(1,length(treated)),rep(0,length(control)))
	
	pvals=hl_max=hl_min=lambda=matrix(0,length(Gamma),length(Delta)) 
	wilx=wilcox.test(x=treated,y=control,exact=F,alternative=alternative)$p.value
	
	qSolve=function(lambda,n,m,r){
		if(lambda!=1){
			a=lambda-1
			b=-((lambda-1)*(m+r)+n)
			c=lambda*r*m
			
			E=c((-b+sqrt(b^2-4*a*c))/(2*a),(-b-sqrt(b^2-4*a*c))/(2*a))
			E=E[which(max(0,r+m-n) <= round(E,digits=10) & round(E,digits=10) <= min(r,m))]
		} else if(lambda==1){
			E=(r*m)/n
		}
		V=1/(1/(E)+1/(r-E)+1/(m-E)+1/(n-r-m+E))
		
		return(list(E=E,V=V))
	}
	
	uSolve=function(E=E,V=V,n=n,m=m,k=k,qsort=qsort,ties=F){
		# recall that qsort is the sorted ranks of the catenated treated and control vectors of responses 
		if(ties==F){
			if(k==0){
				q1=0
				w1=0
				q0=sum(qsort[(k+1):n])/(n-k)
				w0=sum((qsort[(k+1):n]-q0)^2)/(n-k-1)
				mu=sum(qsort)*m/n
				sigma=sqrt((sum((qsort-sum(qsort)/n)^2))*(m*(n-m))/(n*(n-1)))
			} else if(k==n){
				q1=sum(qsort[1:k])/k
				w1=sum((qsort[1:k]-q1)^2)/(k-1)
				q0=0
				w0=0
				mu=sum(qsort)*m/n
				sigma=sqrt((sum((qsort-sum(qsort)/n)^2))*(m*(n-m))/(n*(n-1)))
			} else if(k==1){
				q1=sum(qsort[1:k])/k
				w1=0
				q0=sum(qsort[(k+1):n])/(n-k)
				w0=sum((qsort[(k+1):n]-q0)^2)/(n-k-1)
				mu=E*q1+(m-E)*q0
				sigma=sqrt((w1-w0)*E-(E^2+V)*((w1)/(k)+(w0)/(n-k))+(m*(n-k-m+2*E)*w0)/(n-k)+V*((q1-q0)^2))
			} else if(k==(n-1)){
				q1=sum(qsort[1:k])/k
				w1=sum((qsort[1:k]-q1)^2)/(k-1)
				q0=sum(qsort[(k+1):n])/(n-k)
				w0=0
				mu=E*q1+(m-E)*q0
				sigma=sqrt((w1-w0)*E-(E^2+V)*((w1)/(k)+(w0)/(n-k))+(m*(n-k-m+2*E)*w0)/(n-k)+V*((q1-q0)^2))
			} else if(k>1 & k<(n-1)){
				q1=(sum(qsort[1:k])/k)
				q0=(sum(qsort[(k+1):n])/(n-k))
				w1=sum((qsort[1:k]-q1)^2)/(k-1)
				w0=sum((qsort[(k+1):n]-q0)^2)/(n-k-1)
				mu=E*q1+(m-E)*q0
				sigma=sqrt((w1-w0)*E-(E^2+V)*((w1)/(k)+(w0)/(n-k))+(m*(n-k-m+2*E)*w0)/(n-k)+V*((q1-q0)^2))
			}
		} else if(ties==T){
			mu=(E*n+m*(n-k+1))/2
			sigma=sqrt(((((n+1)*(2*k-n)+2*m*(n-k+1))*E-(n+2)*E^2)/12 + (m*(n-k-m)*(n-k+1))/12 + (V*(n-1)*(n+2/3))/4))
		}
		return(list(mu=mu,sigma=sigma))
	}	
	
	hlSolve=function(R,Z,t=NULL,matched=F,inc=.0001){
		
		m=length(which(Z==1))
		n=length(R)
		invisible(library(stringr))  
		dg=abs(str_length((inc))-3)
		

		#matched == F is rank sum statistic 
		#matched == T is sum sign rank statistic 
		
		if(matched==F & is.null(t)){
			t=m*(n+1)/2
		}
		if(matched==T & is.null(t)){
			t=m*(n/2+1)/4
		}
		
		p_tau=m_tau=0
		np0=np1=nm0=nm1=1
		
		while(sign(np0)==sign(np1) & sign(nm0)==sign(nm1)){
			if(matched==F){
				np0=sum(rank((R-p_tau*Z))[which(Z==1)]) - t
				p_tau=p_tau+inc
				np1=sum(rank((R-p_tau*Z))[which(Z==1)]) - t
				
				nm0=sum(rank((R-m_tau*Z))[which(Z==1)]) - t
				m_tau=m_tau-inc
				nm1=sum(rank((R-m_tau*Z))[which(Z==1)]) - t
				
			} else if(matched==T){
				np0=sum(rank(abs(((R-p_tau*Z)[Z==1])-((R-p_tau*Z)[Z==0])))[which(((R-p_tau*Z)[Z==1])>((R-p_tau*Z)[Z==0]))]) - t
				p_tau=p_tau+inc
				np1=sum(rank(abs(((R-p_tau*Z)[Z==1])-((R-p_tau*Z)[Z==0])))[which(((R-p_tau*Z)[Z==1])>((R-p_tau*Z)[Z==0]))]) - t
				
				nm0=sum(rank(abs(((R-m_tau*Z)[Z==1])-((R-m_tau*Z)[Z==0])))[which(((R-m_tau*Z)[Z==1])>((R-m_tau*Z)[Z==0]))]) - t
				m_tau=m_tau-inc
				nm1=sum(rank(abs(((R-m_tau*Z)[Z==1])-((R-m_tau*Z)[Z==0])))[which(((R-m_tau*Z)[Z==1])>((R-m_tau*Z)[Z==0]))]) - t
			}
		}
		
		if(abs(nm1) + abs(nm0) < abs(np1) + abs(np0)){
			tau=-median(c(-c(m_tau,m_tau+inc),c(p_tau-inc,p_tau)))
		} else if(abs(nm1) + abs(nm0) > abs(np1) + abs(np0)){
			tau=median(c(-c(m_tau,m_tau+inc),c(p_tau-inc,p_tau)))
		} else{
			tau=0
		}
		
		return(tau)	
	}	
	
	for(i in 1:length(Gamma)){ 
		for(j in 1:length(Delta)){ 
			
			lambda[i,j]=(Delta[j]*Gamma[i]+1)/(Delta[j]+Gamma[i])
			
			if(lambda[i,j]==1){
				pvals[i,j]=wilx
				hl_max[i,j]=hlSolve(R=R,Z=Z,t=NULL,matched=F,inc=inc)
			} else {
				mu=sigma=array(NA,n+1)
				for(k in 0:n){
					qS=qSolve(lambda[i,j],n=n,m=m,r=k)
					uS=uSolve(E=qS$E,V=qS$V,n=n,m=m,k=k,qsort=qsort)
					mu[k+1]=uS$mu
					sigma[k+1]=uS$sigma
				} #k-loop
				
				if(lambda[i,j]>1){
					t_max=max(mu)
					pvals[i,j]=1-pnorm(min((t-mu)/sigma))
				} else if(lambda[i,j]<1){
					t_max=min(mu)		
					pvals[i,j]=1-pnorm(max((t-mu)/sigma))
				}
				hl_max[i,j]=hlSolve(R=R,Z=Z,t=t_max,matched=F,inc=inc)
			} # else lambda 
		} # j-loop
	} # i-loop
	
	options(scipen=3)
	colnames(pvals)=colnames(hl_max)=as.character(round(Delta[1:length(Delta)],digits=3))
	rownames(pvals)=rownames(hl_max)=as.character(round(Gamma[1:length(Gamma)],digits=3))
	
	hl_max[1,1]=hl_max[1,1]*sgn
	outs=list()
	outs[[1]]=outs[[2]]=outs[[3]]=lambda=matrix(NA,length(Gamma),length(Delta))
	outs[[1]]=hl_max
	outs[[3]]=outs[[2]]=pvals

	for(j in 1:ncol(outs[[3]])){
		inds1=which(outs[[3]][,j]>.5)
		inds2=which(outs[[3]][,j]<=.5)
		outs[[3]][inds1,j]=2*(1-outs[[3]][inds1,j])		
		outs[[3]][inds2,j]=2*(outs[[3]][inds2,j])
	}

	return(list(HL=round(outs[[1]],digits=digits),pvalues=list(one_tailed=round(outs[[2]],digits=digits),two_tailed=round(outs[[3]],digits=digits)),Gamma=Gamma,Delta=Delta))
} 

#END
